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25  Abstract 

26  In  this  paper,  we  present  a  dynamical  core  for  the  atmospheric  primitive  hydrostatic 

27  equations  using  a  unified  formulation  of  spectral  element  (SE)  and  discontinuous  Galerkin 

28  (DG)  methods  in  the  horizontal  direction  with  a  finite  difference  (FD)  method  in  the  radial 

29  direction.  The  CG  and  DG  horizontal  discretization  employs  high-order  nodal  basis  functions 

30  associated  with  Lagrange  polynomials  based  on  Gauss-Lobatto-Legendre  (GLL)  quadrature 

31  points,  which  define  the  common  machinery.  The  atmospheric  primitive  hydrostatic 

32  equations  are  solved  on  the  cubed-sphere  grid  using  the  flux  form  governing  equations  in  a 

33  three-dimensional  (3D)  Cartesian  space.  By  using  Cartesian  space,  we  can  avoid  the  pole 

34  singularity  problem  due  to  spherical  coordinates  and  this  also  allows  us  to  use  any 

35  quadrilateral-based  grid  naturally.  In  order  to  consider  an  easy  way  for  coupling  the  dynamics 

36  with  existing  physics  packages,  we  use  a  FD  in  the  radial  direction.  The  models  are  verified 

37  by  conducting  conventional  benchmark  test  cases:  the  Rossby-Haurwitz  wavenumber  4, 

38  Jablonowski-Williamson  tests  for  balanced  initial  state  and  baroclinic  instability,  and  Held- 

39  Suarez  tests.  The  results  from  those  tests  demonstrate  that  the  present  dynamical  core  can 

40  produce  numerical  solutions  of  good  quality  comparable  to  other  models.. 

41 
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42  1.  Introduction 


43  Spectral  element  (SE;  here  after  is  referred  to  as  continuous  Galerkin  (CG))  and 

44  discontinuous  Galerkin  (DG)  methods  are  very  attractive  on  many-core  computing  platforms 

45  because  these  methods  decompose  the  physical  domain  into  smaller  pieces  having  a  small 

46  communication  footprint.  CG/DG  methods  are  local  in  nature  and  thus  can  have  a  large  on- 

47  processor  operation  count  (Kelly  and  Giraldo,  2012)  which  is  advantageous  on  large 

48  processor-count  computers.  Also  CG/DG  methods  can  achieve  high-order  accuracy  because 

49  the  polynomial  order  can  be  adjusted  automatically  according  to  the  corresponding  numerical 

50  integration  rule,  that  is,  the  Gaussian  quadrature  (Taylor  et  al.  1997;  Giraldo  2001;  Giraldo  et 

51  al.  2002).  In  addition,  CG/DG  methods  are  geometrically  flexible  in  the  types  of  grids  they 

52  can  use;  this  includes  static  and  adaptive  grids  as  well  as  conforming  and  non-conforming 

53  grids  (Giraldo  et  al.  2002;  Giraldo  and  Rosmond  2004;  Mueller  et  al.  2013). 

54  The  CG  method  is  characterized  by  the  high-order  approximation  combined  with  the 

55  local  decomposition  property  of  the  finite  element  method  (FEM)  and  weak  numerical 

56  dispersion  property  of  the  spectral  method.  The  DG  method,  on  the  other  hand,  is  best 

57  characterized  as  a  combination  of  the  properties  of  the  CG  method  plus  the  local  conservation 

58  properties  of  the  finite  volume  method  (FVM)  (Giraldo  and  Restelli  2008).  The  virtues  of 

59  the  DG  method  are  that  it  is  inherently  conservative  (both  locally  and  globally)  as  in  the  case 

60  of  the  FVM.  However,  the  common  criticism  of  the  DG  method  is  the  stringent  Courant- 

61  Friedrichs-Lewy  (CFL)  stability  constraint  in  explicit  time  schemes.  For  a  DG  method  using 

62  k-th  order  basis  functions,  an  approximate  CFL  limit  estimate  is  l/(2k+l)  (Cockburn  and  Shu 

63  1989).  This,  however,  is  partly  due  to  the  choice  of  the  numerical  flux  which,  for  expediency, 

64  is  chosen  as  a  purely  edge-based  flux  although  other  fluxes  are  also  possible  (e.g.,  Yelash  et 

65  al.  2014);  however  these  more  sophisticated  approaches  come  at  a  price  and  it  is  yet  unclear 
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66  which  strategy  yields  a  faster  wallclock  time  to  solution. 

67  To  date,  successful  applications  of  the  CG  method  in  hydrostatic  atmospheric  modeling 

68  include  the  Community  Atmosphere  Model  -  spectral  element  dynamical  core  (CAM-SE) 

69  (Dennis  et  al.  2012)  and  the  scalable  spectral  element  Eulerian  atmospheric  model  (NSEAM) 

70  (Giraldo  and  Rosmond,  2004,  hereafter  GR04).  In  this  context,  one  of  the  motivations  of  this 

71  study  is  to  construct  a  dynamical  core  using  a  unified  formulation  of  CG  and  DG  methods  as 

72  described  in  Giraldo  and  Restelli  2008  and  Kelly  and  Giraldo  2012  where  nonhydrostatic 

73  atmospheric  models  are  proposed.  Successful  applications  of  the  DG  method  in  hydrostatic 

74  atmospheric  modeling  include  the  work  of  Nair  et  al.  2009;  however,  in  our  paper  we  shall 

75  present  results  for  more  than  one  test  case.  To  our  knowledge,  the  results  for  the  Held-Suarez 

76  test  cases  presented  in  our  paper  are  the  first  such  results  shown  for  a  DG  model.  The 

77  significance  is  that  this  confirms  the  long-tenn  stability  of  the  DG  method  for  hydrostatic 

78  models.  Although  we  could  also  discretize  the  vertical  direction  with  CG  and  DG  methods, 

79  we  choose  a  conservative  flux-form  finite-difference  method  for  discretization  in  the  vertical 

80  direction  which  is  similar  to  the  approach  used  in  both  CAM-SE  and  NSEAM.  This  choice  of 

81  vertical  discretization  provides  an  easy  way  for  coupling  the  dynamics  with  existing  physics 

82  packages. 

83  In  this  paper  we  construct  a  unified  formulation  of  CG  and  DG  for  the  primitive 

84  hydrostatic  equations  in  GR04.  In  order  to  achieve  a  unified  formulation,  the  advective-form 

85  governing  equations  in  GR04  are  recast  in  flux  form.  GR04  provides  a  clue  for  converting  the 

86  advective-form  equation  set  in  3D  Cartesian  space  to  the  flux  fonn  in  their  appendix.  By 

87  using  3D  Cartesian  space,  we  can  be  free  from  the  pole  singularity  problem  in  spherical 

88  coordinates.  Although  a  local  Cartesian  coordinate  system  could  also  be  used  to  overcome 

89  these  problems  (Taylor  et  al.  1997;  Nair  et  al.  2005),  the  use  of  3D  Cartesian  space 

90  everywhere  allows  us  to  treat  the  pole  as  any  other  point.  Therefore  it  permits  general  grids 
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naturally  such  as  icosahedral,  hexahedral,  and  adaptive  unstructured  grids  (it  should  be  noted 
that  general  grids  can  also  be  used  with  the  coordinate  invariant  form  of  the  equations).  In 
this  paper  we  adopt  a  hexahedral  grid  -  the  so  called  cubed-sphere. 

In  brief,  the  objective  of  this  paper  is  to  show  the  feasibility  of  the  hydrostatic  primitive 
equation  models  using  CG/DG  horizontal  discretization  and  the  FD  vertical  discretization  by 
conducting  conventional  benchmark  test  cases.  The  organization  of  the  remainder  of  this 
paper  is  as  follows.  In  the  next  section  we  describe  the  governing  equations  in  3D  Cartesian 
space  with  a  definition  of  the  prognostic  and  diagnostic  variables.  In  Sec.  3  we  explain  the 
horizontal,  vertical,  and  temporal  discretization  methods  including  the  numerical 
approximation  of  the  equations.  In  Sec.  4  we  describe  the  cubed-sphere  grid,  and  in  Sec.  5, 
we  present  the  simulation  results  of  the  test  cases.  Finally,  in  Sec.  6,  we  end  the  paper  with  a 
summary  of  our  findings  and  some  concluding  remarks. 


2.  Governing  Equations 

The  primitive  hydrostatic  equations  of  conservation  form  in  the  3D  Cartesian  space  with 
a  sigma  pressure  vertical  coordinate  a  are  given  as 

37+v'F-s-+VS„  (i) 
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112  respectively  denote  Coriolis  with  the  Lagrange  multiplier  p ,  horizontal,  and  vertical  source 


113  terms,  and 


115  is  the  horizontal  flux  terms  where  /  ,  j  ,  and  k  denote  the  Cartesian  directional  unit 

116  vectors.  The  prognostic  variables  q  are  comprised  of:  1)  the  surface  pressure  n  defined  as 

117  7T  =  Ps~pt,  (5) 

118  where  ps  is  the  true  surface  pressure,  and  pt  is  the  pressure  at  the  top  of  the  atmosphere;  2) 

119  the  flux-form  velocity  components  U  =  (U,V,W)  =  [nu , 7N , ttsn) ,  where  (u,v,w)  are  the 

120  three  Cartesian  velocity  components,  and  3)  the  flux-form  potential  temperature  0  =  nd , 

121  where  0  is  the  potential  temperature.  The  diagnostic  variables  are  1)  the  geopotential  (f) 

122  given  by  the  diagnostic  equation  as 
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(6) 


d  (7  p  —  p  r  - 1 

128  the  a  -coordinate  vertical  velocity  -  where  <j  = - -e  0,1  is  the  definition 

dt  n  L  J 

129  of  the  sigma  pressure  coordinate  with  a  value  of  0  at  the  top  of  the  atmosphere  and  1  at  the 

130  surface.  The  constants  o  and  co  in  Eq.  (4)  are  the  Earth’s  radius  and  angular  velocity, 

131  respectively,  and  //  is  a  Lagrange  multiplier  for  the  fluid  particles  to  remain  on  a  spherical 


132  shell  with  constant  a .  The  momentum  variables  representing  the  atmospheric  motion  over 

133  the  shell  in  the  Cartesian  space  have  three  components  along  the  x,  y,  and  z  axes  in  Cartesian 

134  coordinates,  so  that  the  movement  of  a  particle  on  the  shell  has  three  degrees  of  freedom, 

135  which  can  move  freely  in  R\  To  ensure  that  fluid  particles  remain  on  the  spherical  shell,  it  is 

136  required  that  the  fluid  velocity  remains  perpendicular  to  the  position  vector,  which  yields  a 

137  Lagrange  multiplier  in  the  momentum  equations  (Giraldo  2001;  Giraldo  et  al.  2002;  Giraldo 

138  and  Rosmond,  2004).  It  is  noteworthy  that  among  the  independent  variables  (x ,  y ,z ,(7,t) , 


139  (x ,  y  ,z  )  represent  grid  points  on  the  sphere  which  are  related  to  the  points  in  the  spherical 


140  coordinates  (A,  (p)  given  as 

x  =  o  cos  /Icos  cp, 

141  y  =  o  s in  Acos  (p,  (8) 

z  =  a  sin <p. 

142  Thus  V  is  defined  as 
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143 


(9) 


144  at  constant  a. 

145 

146  3.  Discretization 


147  1)  Discretization  in  the  horizontal  direction 

148  To  describe  the  discretization  of  the  horizontal  operators  by  the  CG/DG  method  we 

149  follow  the  description  given  previously  in  Giraldo  and  Restelli  2008  and  in  Kelly  and  Giraldo 

150  2012.  Let  us  begin  by  rewriting  Eq.  (1)  as  follows 

151  —  +  V  ■  F  =S  (10) 

dt 

152  Next,  let  us  introduce  the  following  vector  spaces 

153  VCG  ={yseH1(Q.)\ysePN(n)}  (11) 

154  and 

155  \/DG  ={yseL2(Q.)\ysePN(n  )}  (12) 

156  where  we  now  seek  solutions  of  Eq.  ( 1)  as  follows: 

157  q  eV  V 

158  where  V  denotes  either  VCG  or  V  DG  .  Next,  we  approximate  the  solution  vector  as  follows 

M 

159  qN(x,y,z,t)  =^l//.(x,y,z)q.(t)  (13) 

/= 1 

160  where,  for  quadrilateral  elements  in  the  horizontal  direction,  M  =  (N  +1)2  with  N 


8 


161 

162 

163 

164 


representing  the  polynomial  order  of  the  basis  function  y/  . 

We  now  introduce  this  expansion  into  our  governing  system  of  equations,  multiply  by  a 
test  function,  and  integrate  by  parts  to  yield 

L  ¥>  +  Jr  ¥"  FdTe-  L  V  V F  )d  =  JQ  '  (14) 

e  C/L  1  e  e 
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where  the  terms  with  Qe  refer  to  volume  integrals  and  the  one  with  Te  is  a  boundary 

integral  which  accounts  for  both  internal  faces  (neighboring  elements  share  faces)  as  well  as 
boundary  faces  (elements  on  boundaries  do  not  share  faces  with  other  elements).  In  matrix- 
vector  form,  this  equation  can  be  written  as 


where 


k  ^ + H;  )r  f>»  )  -  )r 


M'l 

K;e  =  lr  ¥¥Jdd  V?’ 

dt.  =  [  V v/.v/.d  Q  . 

'■>  Jn  r'rJ  e 


(15) 


(16) 


These  matrices  represent:  the  mass,  flux,  and  differentiation  matrices,  respectively. 

For  the  DG  method,  the  matrix-vector  fonn  given  above  is  sufficient  as  long  as  we  define 
the  numerical  flux  ,  e.g.,  as  follows 


F  *(Q, 


(qLN)+F(qRN)~n\Amax\(qRN  -qLN) 


(17) 


where  the  superscripts  L  and  R  refer  to  the  left  and  right  elements  (arbitrarily  decided)  of  the 
face  Te  and  A  is  the  maximum  eigenvalue  of  the  Jacobian  matrix  of  the  governing 

partial  differential  equations.  Here  we  use  the  Rusanov  scheme  for  the  numerical  flux 
because  of  its  simplicity  although  any  other  Riemann  solver  could  be  used.  For  the  CG 
method,  the  matrix-vector  form  given  above  is  also  used  except  that  the  tenn  of  the  flux 
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181  matrix  vanishes  on  the  sphere  and  we  then  use  the  direct  stiffness  summation  (DSS)  operation 

182  which  gathers  the  element-wise  solution  to  a  global  grid  point  solution  and  then  scatters  it 

183  back  to  the  element-wise  space.  This  is  done  to  ensure  that  the  solution  is  C°  across  all 

184  element  faces. 

185 

186  2)  Discretization  in  the  vertical  direction 

187  We  use  the  FD  method  similarly  to  other  global  models  to  gain  an  easy  way  for  coupling 

188  the  dynamics  with  existing  physics  packages,  although  we  could  also  discretize  the  vertical 

189  operators  with  the  CG/DG  methods  (as  done  in  Kelly  and  Giraldo  2012;  Giraldo  et  al.  2013). 

190  Also  by  using  the  FD,  we  can  keep  the  model  as  similar  as  possible  to  the  NSEAM  model 

191  (GR04)  so  that  we  directly  discern  differences  from  the  discrete  horizontal  operators.  Using  a 

192  Lorenz  staggering,  the  variables  U  ,  V  ,  W  ,  0  ,  and  (p  are  at  layer  mid  points  denoted  by 

193  k  —  Nlev  where  Nlev  is  the  total  number  of  layers,  while  the  variable  P  and 

1 

194  are  at  layer  interface  points  denoted  by  k  +  —,  k  -  0,X...,Nlev  . 

d/r 

195  We  begin  the  vertical  discretization  by  the  evaluating  —  which  is  given  by 

3t 

196  integrating  the  first  row  of  Eq.  (1)  (i.e.,  the  continuity  equation)  from  the  surface 

197  (<r  =  <j„,  =1)  to  the  top  (<r  =  <t.  =  0)  with  no-flux  boundaries  at  the  top  and 

v  bottom  Nlev+V  2  7  r  v  top  1/2  7  r 

198  bottom  levels  of  the  atmosphere  (i.e.  dkp  =  <£fcottom  =0).  Thus, 

A  jj.  Nlev 

199  - =  IVT1  <7t,  (18) 

ot  k=1 

200  where  k  is  the  number  of  vertical  levels  to  be  integrated  across  and  Acr;  =  <J/+V2  -  tJ,_V2  is 

201  the  thickness  of  the  layer.  Then  the  vertical  velocity  at  each  vertical  level  is  obtained  by 

202  integrating  the  continuity  equation  from  the  top  of  the  atmosphere  to  the  material  surface  as 
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203  follows 


204 


207 


212 


r)7T  k  ■ 

((77r)«+i/2  =  -  — ok+V2  +  M-'X  D  ■  . 


(19) 


/= i 


205  The  vertical  advection  term  — — -  in  the  vertical  source  term  Sv  is  computed  using  the 

206  third-order  upwind  biased  discretization  in  Hundsdorfer  et  al.  (1995)  which  is  given  as 


da 


_  4 -2  -  84  -1  +  84  +1  +2  !  3  jg  n(  ^4 -2  ~  44  -1  +  64  ~  44  +1  +4 +2 


12A<7 


12A<t 


(20) 


208  where  /  denotes  the  flux  (<&/).  It  is  noted  that  the  upwind-biased  schemes  are  inherently 

209  diffusive.  Following  GR04,  the  hydrostatic  equation,  Eq.  (6),  is  evaluated  as  follows 


210  fik +1  C  p®k^k+l]2  +  C  p  +1^*-  +1  Pfr+1/2)’ 

211  where  the  Exner  function  at  layer  interfaces  and  midpoints  is  given  by 


k  +1/2 


(  Y 
Pk+y  2 

P  0 


(21) 


(22) 


213  and 


214 


215  respectively. 

216 


217 

218 

219 

220 

221 


1  1 


fpK+1 


K  +  1Po 


p ' 


k  +3/2  r  k  -1/2 
yP k+yi  —  P k-y 2  j 


(23) 


3)  Discretization  in  time 

For  integrating  the  equations,  we  adopt  a  third-order  strong  stability  preserving  explicit 
Runge-Kutta  (SSP-RK)  scheme  (Cockbum  and  Shu  1998;  Nair  et  al.  2005).  The  3rd  order 
SSP-RK  scheme  is  introduced  into  our  governing  equations  in  the  form  of 


dt 


(24) 
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222  and  is  given  as  follows: 


q[1)  =qn  +  AtR(qn) 

223  q <2)  =  —q  n  +  —q  (1)  +  —  At  R  (q(1))  (25) 

4  4  4 

<f*‘  =  V  +jq'"+l  AtR(q'!l), 

3  3  3 

224  where  the  superscripts  n  and  n  + 1  denote  time  levels  t  and  t  +  At ,  respectively.  While 

225  for  smooth  problems  the  SSP-RK  scheme  does  not  generate  spurious  oscillations  so  that  are 

226  widely  used  for  DG  methods,  for  problems  with  strong  shocks  or  discontinuities,  oscillations 

227  can  lead  to  nonlinear  instabilities  (Cockburn  and  Shu  1998).  Since  an  SSP-RK  time- 

228  integration  scheme  cannot  control  such  undesirable  effects,  a  Boyd-Vandeven  spatial  filter  is 

229  applied  after  the  time  integration,  which  is  described  in  GR04.  Neither  viscosity  nor  slope 

230  limiter  are  used  in  all  simulations. 

231 

232  4.  Cubed-sphere  Grid 

233  The  cubed-sphere  grids  are  composed  of  the  six  patches  obtained  by  the  gnomonic 

234  projection  of  the  faces  of  the  hexahedron  which  are  subdivided  into  (nHxnH)  quadrilateral 

235  elements  where  nH  is  the  number  of  quadrilateral  elements  in  each  direction  (GR04).  Inside 

236  each  element  we  build  (N  + 1)  Gauss-Lobatto-Legendre  (GLL)  quadrature  points,  where  N 

237  indicate  the  polynomial  order  of  the  basis  function  y/  .  Therefore  the  total  number  of  grid 

238  points  N  is  given  as 

239  Np  =6(nHA/)2  +2,  (26) 

240  and  the  number  of  elements  Ne  comprising  the  sphere  is 

241  Ne=6(nHf.  (27) 
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242  We  now  introduce  the  square  region  on  the  gnomonic  space  (gG ,  rjG  j  = 


246 

247 

248 

249 

250 

251 


253 


254 


k  n 
,3 — 

4  4 


in 


243  each  of  the  six  faces  to  describe  the  relation  to  spherical  coordinates  (A,  rpj .  The  gnomonic 


244  space  (4,/7g)  = 

245  (4'%  )  via 


k  n 
— 

4  4 


is  mapped  to  the  corresponding  spherical  coordinates 


4  =4 


f 


cp  =  a  res  in 


tan  T) 


(28) 


(29) 


/ l+tan24  +tan2^G  ^ 

and  then  we  construct  the  cubed-sphere  grid  by  rotating  this  face  to  the  six  faces  of  the 
hexahedron  by 


A  =  A  +  a  rc  ta  n 


cos  (pG  sin  Ag 


c os  cpG  c os  Ag  c os  (pc  —  sin  (pG  s in  (pc 


(p  =  a  res  in  (s  in  cpG  cos  (pc  +cos  (pG  cos  AG  sin  (pc 


252  with  the  centroids,  (Ac,  (pc )  =  [c  -l]— ,0  for  c  K,4  ,  (yAs,cps)  = 


°-i 

v  / 


(30) 

(31) 

,  and 


{\’<Pe)  = 


1 

°'-I 

V  J 


The  resolution  of  the  cubed-sphere  grid  H  is  detennined  by  n  (the  number  of 


255  quadrilateral  elements  in  each  direction  contained  in  each  of  the  six  faces  of  the  cube)  and  N 

256  (the  polynomial  order  of  the  elements),  where  we  use  H  =nHN  as  the  convention  to  define 

257  the  grid  resolution.  Fig.  1  show  examples  of  the  grids  with  H  =  3  (n  =  3  and  N  =1), 


258  H  =15  (n  =3  and  N  =  5 ),  and  H  =35  (n  =5  and  N  =  7 ). 
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259 


260  5.  Simulation  results  with  Benchmark  Tests 

261  We  consider  the  following  test  cases:  1)  3D  Rossby-Haurwitz  wavenumber  4, 

262  2)  Jablonowski-Williamson  balanced  initial  state  test,  3)  baroclinic  instability  test,  and  4) 

263  Held-Suarez  test.  Because  all  of  the  test  cases  except  2)  the  Jablonowski-Williamson 

264  balanced  initial  state  test  do  not  have  analytical  solutions,  we  compare  our  results  to  the 

265  results  of  other  published  papers  and  evaluate  the  results  qualitatively.  We  now  discuss  the 

266  results  of  the  four  test  cases. 

267 

268  1)  3D  Rossby-Haurwitz  wavenumber  4 

269  We  conduct  the  Rossby-Haurwitz  (RH)  wave  test  case  which  is  a  3D  extension  of  the 

270  2D  shallow  water  RH  wave  discussed  in  Williamson  et  al.  (1992).  The  main  differences 

271  compared  to  the  2D  shallow  water  formulation  include  the  introduction  of  a  temperature  field 

272  and  the  derivation  of  the  surface  pressure,  which  is  discussed  in  GR04  and  Jablonowski  et  al. 

273  (2008).  The  Rossby-Haurwitz  wave  approximately  preserves  its  shape  even  in  nonlinear 

274  shallow  water  and  primitive  equation  models,  which  has  a  sufficiently  simple  enough  pattern 

275  to  allow  one  to  judge  if  the  simulation  was  successful.  We  initialize  the  model  following 

276  Jablonowski  et  al.  (2008). 

277  Snapshots  of  the  output  data  for  the  CG  and  DG  models  for  day  15  are  presented  in  Figs. 

278  2  and  3,  respectively.  The  figures  show  the  850  hPa  zonal  wind,  meridional  wind,  and 

279  temperature  as  well  as  the  surface  pressure.  These  model  results  were  computed  at  the 

280  resolution  of  H  =  64  (nH  =  8  and  N  =  8  )  with  26  vertical  levels  (Nlev=26).  The  results 

281  of  the  CG  and  DG  simulations  are  virtually  indistinguishable;  in  addition,  the  accuracy  results 

282  of  both  simulations  are  almost  identical  to  the  results  obtained  with  the  CAM3.5.41  version 
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of  the  NCAR  Finite  Volume  (FV)  dynamical  core  at  the  resolution  1°  by  1°  with  26  hybrid 
levels,  as  described  in  Jablonowski  et  al.  (2008).  Although  we  have  used  a  relatively  low 
resolution  of  H64  which  is  comparable  to  T63  of  a  spectral  model,  the  results  are  strikingly 
similar  to  the  solutions  with  the  l°xl°  NCAR  CAM-FV  core,  both  in  phase  and  amplitude. 


2)  Jablonowski-Williamson  balanced  initial  state  test 

In  order  to  estimate  the  accuracy  and  stability  of  the  dynamical  core,  we  conduct  the 
Jablonowski-Williamson  balanced  initial  state  test  introduced  by  Jablonowski  and 
Williamson  (2006).  We  initialize  the  model  following  Jablonowski  and  Williamson  (2006a 
and  b).  Using  the  balanced  initial  fields,  the  simulation  results  should  maintain  the  initial  state 
perfectly  for  a  sufficient  amount  of  time.  Since  the  initial  state  of  this  test  is  the  true  solution, 
we  can  compute  error  norms.  We  evaluate  the  error  by  using  the  relative  L2  error  defined  by 


simulation 


f  (q  ,  -q  ...  fd  Q. 

JQ.\  exact  'simulation  ) 


L 


dQ. 


where  q  . represents  the  computed  state  variables  and  q  ,  the  exact  (i.e.,  initial 

'  simulation  r  r  1  exact  v  7 

condition)  values. 

Figure  4  shows  the  normalized  surface  pressure  L2  error  norms  for  the  CG  and  DG 
simulations  with  A7  =  1  28  (n H  =16  and  N  =  8 )  horizontal  resolution  and  26  vertical 

levels  (Nlev=26).  The  L2  error  norms  of  the  two  simulations  are  visually  identical,  in  which 
the  error  oscillates  but  remains  bounded.  These  results  (including  the  value  of  the  L2  error) 
compare  well  against  those  of  the  NSEAM  model  presented  in  GR04.  The  bounded  error 
confirms  that  the  initial  balanced  state  is  properly  maintained.  In  practice  though,  the  initial 
state  degrades  over  time.  After  20-days,  the  zonal  wind  fields  for  the  CG  and  DG  simulations 
show  a  somewhat  distorted  distribution  with  an  increasing  zonally  asymmetric  pattern  (Fig. 
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5).  Initially  the  maximum  of  the  zonal  winds  at  the  lowest  level  are  about  9.4  m/s  in  mid¬ 
latitude,  but  after  20-days  the  maximum  difference  of  the  zonal  wind  is  up  to  about  0.02  m/s 
showing  the  zonal  asymmetry.  Although  the  error  distribution  is  different  between  the  CG 
and  DG  simulations  in  detail,  these  have  a  wavenumber  4  structure  which  arise  from  the 
cubed-sphere  grid.  The  wavenumber  4  signals  grow  over  time  and  lead  eventually  to  a 
breakdown  of  the  balanced  state.  However,  higher  resolutions  delay  the  growth  of  the  signals 
as  the  truncation  error  associated  with  the  spatial  discretization  decreases.  Actually,  at 
H  =  192  (nH  =16  and  N  =12)  horizontal  resolution  this  error  virtually  disappears  for 
20-day  simulations  (Fig.  6). 


3)  Jablonowski-Williamson  baroclinic  instability  test 

The  baroclinic  instability  test  case  starts  from  the  balanced  initial  fields,  which  is 
described  above,  with  a  perturbation  in  the  initial  zonal  velocity.  The  baroclinic  wave  is 
induced  by  the  small  perturbation  in  the  initial  zonal  wind.  Here  a  Gaussian  profile  is  used  for 


the  zonal  wind  perturbation,  which  is  centered  at  (/l  ,tpc)  = 
location^20°E,40°N  j.  This  perturbation  is  given  by 


(  K  2k  ^ 

9'  9 

v  j 


pointing  to  the 
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U  perturbation^' (P'(J)=eX  P 


r 

A 

2 

r 

R 

V 

J 

323  where 

324  r  =  a  a rccos  [^s in (pc  sin  (p  +  cos  (pc  cos  ^>cos  (/l  -  A  )  , 

325  and  R  =o/10  is  the  perturbation  radius  (Jablonowski  and  Williamson 2006a  and  b). 

326  Since  the  baroclinic  wave  test  case  does  not  have  an  analytic  solution,  we  compare  our 

327  results  to  the  solutions  from  Jablonowski  and  Williamson  (2006a)  and  the  NSEAM  model  in 
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328  GR04.  We  show  the  surface  pressure,  850  hPa  temperature,  and  850  hPa  relative  vorticity  at 

329  day  9  for  the  CG  and  DG  simulations  with  the  resolution  of  H  =  80  (n  =16  and 

330  N  =  5  )  and  26  vertical  levels  (Nlev=26)  in  Fig.  7  which  can  be  compared  with  the  solutions 

331  of  the  National  Center  for  Atmospheric  Research’s  Community  Atmosphere  Model  version  3 

332  (NCAR  CAM3)  Eulerian  dynamical  core  at  T85  resolution  and  finite  volume  core  at  1°  by 

333  1.25°  from  Jablonowski  and  Williamson  (2006a).  The  CG  and  DG  simulations  in  Fig.  7  are 

334  visually  very  similar  to  those  reported  in  Jablonowski  and  Williamson  with  regard  to  the 

335  structure  in  the  fields  and  the  extrema  for  the  surface  pressure;  in  addition,  the  CG  and  DG 

336  results  are  almost  identical  to  each  other.  Differences,  however,  can  only  be  seen  in  the 

337  relative  vorticity  field  at  very  small  scales.  In  the  CG  simulation,  the  small-scale  vorticity  in 

338  the  vicinity  of  the  hook  is  depicted,  and  the  maximum  strength  of  the  relative  vorticity  is 

339  larger  than  that  of  the  DG  simulation,  which  can  be  also  seen  in  the  results  of  a  relatively 

340  higher  resolution  shown  in  Fig.  8.  Figure  8  shows  the  same  fields  at  day  9  as  in  Fig.  7  but  for 

341  the  higher  resolution  of  H  =160  ( nH  =32  and  N  =  5  )  and  26  vertical  levels  (Nlev=26). 

342  In  comparison  with  the  results  of  the  lower  resolution  of  H  =  80  (n  =16  and  N  =  5  ),  it 

343  can  be  clearly  seen  that  the  numerical  solutions  of  the  two  different  resolutions  are  well 

344  converged  in  terms  of  the  strength  and  structure  in  the  surface  pressure,  temperature,  and 

345  vorticity  fields.  It  is  noted  that  the  vorticity  fields  in  the  higher  resolution  are  characterized  by 

346  the  smallest  scale  in  the  vicinity  of  the  hook,  which  is  the  same  as  in  the  lower  resolution, 

347  which  imply  that  the  DG  simulation  is  more  diffusive  than  the  CG  simulation.  It  suggests  that 

348  the  diffusive  property  of  the  DG  simulation  is  induced  by  the  Rusanov  numerical  flux  used  in 

349  this  study,  because  the  only  difference  between  the  CG  and  DG  fonnulations  is  the  numerical 

350  flux  and  the  fact  that  the  DG  solutions  are  allowed  to  contain  jumps  across  element  edges. 

351  However,  this  difference  in  the  results  suggests  that  it  is  the  dissipation  of  the  numerical  flux 
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352  that  is  mainly  responsible  for  the  differences  in  the  two  simulations. 

353  In  general,  the  baroclinic  wave  grows  observably  around  day  4.  At  day  7  the  baroclinic 

354  wave  evolves  rapidly  and  by  day  9  the  wave  train  has  intensified  significantly  (Jablonowski 

355  and  Williamson  2006a).  In  order  to  examine  the  growth  of  the  perturbation,  an  evolution  of 

356  the  minimum  surface  pressure  is  shown  in  Fig.  9  which  we  now  compare  with  the  results  in  G 

357  R04.  The  results  of  the  CG  and  DG  simulations  with  different  resolutions  are  almost  in 

358  agreement  until  day  10,  at  which  point  the  simulations  begin  to  show  slight  deviations  from 

359  each  other.  The  DG  simulation  with  the  lower  resolution  tends  to  simulate  somewhat  weak 

360  deepening.  During  the  period  between  day  10  and  11  when  wave  breaking  has  set  in,  the 

361  remarkable  weak  deepening  is  shown  in  the  DG  simulation  at  the  lower  resolution.  At  day  14, 

362  the  difference  of  the  minimum  surface  pressure  between  the  DG  simulation  at  the  lower 

363  resolution  and  the  three  other  simulations  is  about  2  hPa. 

364 

365  4)  Held-Suarez  test 

366  In  order  to  estimate  the  capabilities  of  the  model  in  simulating  a  realistic  climate 

367  circulation  without  complex  parameterizations,  we  conduct  the  Held-Suarez  test.  The  Held- 

368  Suarez  test  ensures  that  a  dynamical  core  produces  a  realistic  zonal  and  time  mean  climate 

369  and  synoptic  eddies  by  using  a  simple  Newtonian  relaxation  of  the  temperature  field  and  a 

370  Rayleigh  damping  of  low-level  winds  representing  boundary-layer  friction  (Held  and  Suarez 

371  1994).  The  Newtonian  relaxation  of  the  temperature  is  added  as  the  diabatic  forcing  term  to 

372  the  thermodynamic  equation,  the  fifth  row  of  Eq.  (1),  and  the  Rayleigh  damping  is  imposed 

373  as  dissipation  tenn  in  the  momentum  equation,  the  second  to  fourth  rows  of  Eq.  (1).  The 

374  detailed  specifications  are  adapted  from  Held  and  Suarez  (1994).  For  this  test  we  use  a 

375  relatively  low  resolution  of  H  =  40  {  nH  =  8  and  N  =  5  )  with  25  vertical  levels 

(Nlev=25)  because  this  test  case  requires  a  relatively  long  model  time  simulation  for  1200 
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376 


377  days.  In  this  paper,  the  integrations  start  from  a  stably  stratified  state  at  rest  atmosphere,  in 

378  which  the  lapse  rate  of  temperature  is  6.5  K/m  and  the  surface  temperature  is  288  K.  We  use 

379  the  simulation  results  from  day  200  to  day  1200  integrations  sampled  every  10-days. 

380  Fig.  10  shows  the  time  mean  zonally  averaged  zonal  wind  and  temperature  for  both  the 

381  CG  and  DG  simulations  which  can  be  easily  compared  to  the  results  of  other  published 

382  papers.  In  comparison  with  the  results  of  the  spectral  transform  model  in  Held  and  Suarez 

383  (1994),  both  the  CG  and  DG  simulations  show  reasonable  and  comparable  distributions, 

384  where  the  midlatitude  jets  at  the  upper  troposphere  near  250  hPa  and  the  equatorial  easterly 

385  flow  in  the  lower  and  upper  atmosphere  are  clearly  visible  in  each  hemisphere.  Also 

386  temperature  stratification  is  maintained  realistically.  The  simulation  results  are  comparable  to 

387  that  of  GR04.  There  exist,  however,  differences  between  the  results  of  the  CG  and  DG 

388  simulations  mainly  in  the  strength  of  the  westerly  flow  and  the  temperature  structure  in  the 

389  upper  atmosphere.  DG  simulates  broader  upper-level  jet  streams  than  CG  that  strengthen  with 

390  altitude.  Also  in  the  temperature  field,  the  DG  simulation  shows  wanner  air  in  the  equatorial 

391  upper  atmosphere.  The  difference  is  shown  clearly  in  Fig.  11  where  we  plot  the  time  mean 

392  zonally  averaged  eddy  heat  flux  of  the  CG  and  DG  simulations.  There  are  two  maxima  at 

393  mid-latitude  in  the  lower  and  upper  atmosphere  indicating  transportations  of  heat  in  the 

394  poleward  direction,  of  which  the  distributions  in  the  CG  and  DG  simulations  are  in  good 

395  agreement  with  previous  studies,  for  example,  Held  and  Suarez  (1994),  Lin  (2004)  and  Wan 

396  et  al.  (2008).  However,  in  comparison  of  the  strength  and  horizontal  gradient  of  the  eddy  heat 

397  flux  between  both  simulations,  CG  simulates  a  stronger  eddy  motion  than  DG. 

398 

399  6.  Summary  and  Conclusions 

400  We  have  proposed  a  hydrostatic  dynamical  solver  using  both  the  continuous  Galerkin 
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401  (CG)  and  discontinuous  Galerkin  (DG)  methods.  It  is  solved  on  a  cubed-sphere  grid  in  3D 

402  Cartesian  coordinates  although  in  principle  any  quadrilateral-based  grid  could  be  used.  The 

403  CG  and  DG  horizontal  discretization  employs  a  high-order  nodal  (Lagrange)  basis  function 

404  based  on  quadrilateral  elements  and  GLL  quadrature  points  which  compose  the  common 

405  machinery.  However,  the  DG  method  use  fluxes  along  the  boundaries  of  the  elements  which 

406  are  approximated  by  the  Rusanov  method.  In  the  vertical  direction,  a  conservative  flux-form 

407  finite-difference  method  is  employed  for  coupling  the  dynamics  with  existing  physics 

408  packages  easily;  we  hope  to  report  progresses  on  this  specific  topic  in  the  future.  A  third- 

409  order  strong  stability  preserving  Runge-Kutta  scheme  was  used  for  time  integration  although 

410  other  time-integrators  (including  semi-implicit  methods)  could  also  be  used. 

411  In  this  paper,  we  show  simulations  of  the  model  using  four  baroclinic  test  cases 

412  including:  the  Rossby-Haurwitz  wave,  balanced  initial  state,  baroclinic  instability,  and  Held- 

413  Suarez  test  cases.  All  cases,  except  for  the  Jablonowski-Williamson  balanced  initial  state  test 

414  case,  do  not  have  analytic  solutions.  Therefore,  we  compare  our  results  to  the  results  of  test 

415  cases  run  by  a  vast  community.  Through  our  comparison  of  the  CG  and  DG  simulations,  we 

416  show  that  for  the  baroclinic  instability  test  and  Held-Suarez  test  cases,  the  DG  simulation 

417  tends  to  simulate  somewhat  weaker  small-scale  features,  such  as  the  minimum  surface 

418  pressure  perturbation  and  eddy  heat  flux,  than  the  CG  method.  This  could  be  due  to  the 

419  intrinsic  diffusion  of  the  Rusanov  numerical  flux  scheme  used  for  the  horizontal 

420  discretization  of  the  DG  method,  which  is  the  only  difference  between  the  CG  and  DG 

421  formulations.  One  of  the  valuable  contributions  of  this  model  is  that  we  can  use  it  to  study  the 

422  effects  of  using  different  horizontal  discretizations  since  we  use  the  exact  same  model  with 

423  the  same  finite  difference  method  in  the  vertical  and  time-integration  methods  but  use  either 

424  CG  or  DG  in  the  horizontal.  The  discrete  operators  in  the  horizontal  use  the  exact  same 

425  numerical  machinery  and  so  the  results  shown  here  isolate  the  differences  offered  by  the  CG 
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426  and  DG  methods.  However,  for  the  other  two  test  cases  (Rossby-Haurwitz  wave  and  balanced 

427  initial  state  tests),  the  results  of  the  CG  and  DG  simulations  are  virtually  indistinguishable. 

428  Furthermore,  the  numerical  results  obtained  for  all  four  test  cases  show  that  the  present 

429  dynamical  core  can  produce  numerical  solutions  of  good  quality  comparable  to  other  models. 

430  The  results  confirm  that  the  CG  and  DG  methods  combined  with  the  finite  difference  method 

431  in  the  vertical  direction  offer  a  viable  strategy  for  atmospheric  modeling.  To  our  knowledge, 

432  we  present  the  first  results  for  a  DG  model  for  long-time  simulations  represented  by  the  Held- 

433  Suarez  test  case.  The  importance  of  this  result  is  that  this  confirms  the  stability  of  the  DG 

434  method  for  long-time  simulations  in  hydrostatic  atmospheric  dynamics.  In  order  to  make  the 

435  model  efficient  and  competitive  with  operational  models,  we  need  a  semi-implicit  time 

436  integration  method  which,  although  requires  some  additional  machinery  to  be  added,  does  not 

437  pose  any  theoretical  barriers  since  such  algorithms  have  already  been  designed  by  one  of  the 

438  authors  in  previous  papers  (Giraldo  2005,  Giraldo  et  al.  2013). 

439 
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534  Figure  Captions 

535  FIG.  1.  The  cubed-sphere  grid  for  (a)  the  H  =  3  (nH  =3  and  N  - 1),  (b)  the 

536  H  =15  (nH  =3  and  N  =  5),  and  (c)  the  H  =35  (nH  =5  and  N  =  7)  horizontal 

537  resolutions. 

538 

539  FIG.  2.  Numerical  results  for  the  CG  simulation  on  the  resolution  of  the  H  =64 

540  ( n  =8  and  N  =  8  )  with  26  vertical  levels:  Top  row:  850  hPa  zonal  wind  and  meridional 

541  wind,  bottom  row:  surface  pressure  and  850  hPa  temperature. 

542 

543  FIG.  3.  As  in  Fig.  2  but  for  the  DG  simulation. 

544 

545  FIG.  4.  L2  error  norm  of  surface  pressure  in  Pa  for  the  CG  and  DG  simulations  at  the 

546  /7  =  1  28  (nH  =16  and  N  =  8  )  horizontal  resolution  and  26  vertical  levels. 

547 

548  FIG.  5.  Distribution  of  zonal  wind  difference  at  the  lowest  model  level  between  day  20 

549  and  day  0  for  the  (top)  CG  and  (bottom)  DG  simulations  at  the  H  =  128  (nH  =16  and 

550  N  =  8  )  horizontal  resolution  and  26  vertical  levels. 

551 

552  FIG.  6.  As  in  Fig.  5  but  for  the  H  =192  (nH  =16  and  N  =12)  horizontal 

553  resolution. 

554 

555  FIG.  7.  Baroclinic  wave  at  day  9  with  the  (left)  CG  and  (right)  DG  simulations  with  the 

556  resolution  of  the  H  =  80  ( nH  =16  and  N  =  5  )  horizontal  resolution  and  26  vertical 
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557  levels:  (upper  row)  surface  pressure,  (middle  row)  850  hPa  temperature,  and  (bottom  row) 

558  850  hPa  relative  vorticity  at  days  (left)  7  and  (right)  9. 

559 

560  FIG.  8.  As  in  Fig.  7  but  for  the  H  =160  ( nH  =32  and  N  -  5). 

561 

562  FIG.  9.  The  minimum  surface  pressure  (hPa)  as  a  function  of  days  for  the  CG  and  DG 

563  simulations  with  the  lower  resolution  of  the  H  =  80  (n  =16  and  N  =  5  )  and  the  higher 

564  resolution  of  the  H  =160  (n  H  =32  and  N  -  5). 

565 

566  FIG.  10.  The  (left)  mean  zonally  averaged  zonal  velocity  (m/s)  and  (right)  mean  zonally 

567  averaged  temperature  (K)  for  the  (upper  row)  CG  and  (bottom  row)  DG  simulations  with  the 

568  resolution  of  the  H  =  40  (n H  =8  and  N  =  5  )  and  25  vertical  levels  (Nlev=25).  These 

569  are  calculated  over  the  last  1000  days  of  a  1200-day  integration. 

570 

571  FIG.  11.  The  mean  zonally  averaged  eddy  heat  flux  for  the  (left)  CG  and  (right)  DG 

572  simulation  with  the  resolution  of  the  H  =40  (n  =8  and  N  =  5). 

573 
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576 

577 

578 


FIG.  1.  The  cubed-sphere  grid  for  (a)  the  H  =  3  (nH  =3  and  N  =1),  (b)  the 
H  =15  (nH  =  3  and  /V  =  5),  and  (c)  the  H  =35  (nH  =5  and  N  =  7)  horizontal 
resolutions. 
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Snapshots  of  the  Rossby-Haurwitz  wave  at  day  15 
simulated  with  KIAPS-H_CG.Np08xNe08L26 


zonal  wind  at  850  hPa  m/s 
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FIG.  2.  Numerical  results  for  the  CG  simulation  on  the  resolution  of  the  H  =64 
(nH  =8  and  N  -  8  )  with  26  vertical  levels:  Top  row:  850  hPa  zonal  wind  and  meridional 
wind,  bottom  row:  surface  pressure  and  850  hPa  temperature. 
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Snapshots  of  the  Rossby-Haurwitz  wave  at  day  1 5 
simulated  with  KIAPS-H_DG.Np08xNe08L26 
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FIG.  3.  As  in  Fig.  2  but  for  the  DG  simulation. 
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589  FIG.  4.  L2  error  norm  of  surface  pressure  in  Pa  for  the  CG  and  DG  simulations  at  the 


590  H-  128  (nH  =16  and  N  -  8  )  horizontal  resolution  and  26  vertical  levels. 
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593  FIG.  5.  Distribution  of  zonal  wind  difference  at  the  lowest  model  level  between  day  20 

594  and  day  0  for  the  (top)  CG  and  (bottom)  DG  simulations  at  the  H  -  1  28  (nH  =16  and 

595  N  -  8  )  horizontal  resolution  and  26  vertical  levels. 
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FIG.  6.  As  in  Fig.  5  but  for  the  H  =192  (n  =16  and  N  =12)  horizontal 


resolution. 
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602  FIG.  7.  Baroclinic  wave  at  day  9  with  the  (left)  CG  and  (right)  DG  simulations  with  the 

603  resolution  of  the  H  =  80  ( nH  =16  and  N  =  5  )  horizontal  resolution  and  26  vertical 

604  levels:  (upper  row)  surface  pressure,  (middle  row)  850  hPa  temperature,  and  (bottom  row) 

605  850  hPa  relative  vorticity  at  days  (left)  7  and  (right)  9. 
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FIG.  9.  The  minimum  surface  pressure  (hPa)  as  a  function  of  days  for  the  CG  and  DG 
simulations  with  the  lower  resolution  of  the  H  =  80  ( n  =16  and  N  =  5  )  and  the  higher 
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FIG.  10.  The  (left)  mean  zonally  averaged  zonal  velocity  (m/s)  and  (right)  mean  zonally 
averaged  temperature  (K)  for  the  (upper  row)  CG  and  (bottom  row)  DG  simulations  with  the 
resolution  of  the  H  =  40  (n  =8  and  N  =  5)  and  25  vertical  levels  (Nlev=25).  These 
are  calculated  over  the  last  1000  days  of  a  1200-day  integration. 
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